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Unlike macroscopic engines, the molecular machinery of living cells is strongly affected by fluctua¬ 
tions. Stochastic Thermodynamics uses Markovian jump processes to model the random transitions 
between the chemical and configurational states of these biological macromolecules. A recently 
developed theoretical framework [Wachtel, Vollmer, Altaner: “Fluctuating Currents in Stochastic 
Thermodynamics I. Gauge Invariance of Asymptotic Statistics”] provides a simple algorithm for the 
determination of macroscopic currents and correlation integrals of arbitrary fluctuating currents. 

Here, we use it to discuss energy conversion and nonequilibrium response in different models for 
the molecular motor kinesin. Methodologically, our results demonstrate the effectiveness of the 
algorithm in dealing with parameter-dependent stochastic models. For the concrete biophysical 
problem our results reveal two interesting features in experimentally accessible parameter regions: 

The validity of a non-equilibrium Green-Kubo relation at mechanical stalling as well as a negative 
differential mobility for superstalling forces. 
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I. INTRODUCTION 


Understanding the complex biochemical processes 
which are responsible for cellular metabolism is one of 
the key questions in modern biophysics. The quanti¬ 
tative analysis of so-called molecular motors, which are 
the small machines transforming different forms of energy 
into one another, is at the center of these efforts m- 
In recent years scientists developed techniques that allow 
the systematic observation and manipulation of these bi¬ 
ological macromolecules [4]. Under in vivo conditions, 
(electro-)chemical gradients in the cell maintain these 
systems out of equilibrium. From a thermodynamic per¬ 
spective, one is interested in the currents of heat, matter 
and energy that flow through a molecular motor, because 
they allow, for instance, the definition of its efficiency. 

In analogy to macroscopic engines, molecular motors 
are described by thermodynamic cycles in a space of bio¬ 
chemical and configurational states. In contrast, the en¬ 
ergy scales involved in biochemical energy conversion are 
only a couple of times larger than the thermal energy. 
Consequently, thermal fluctuations cannot be neglected, 
and their dynamics must be modeled as a stochastic pro¬ 
cess that reproduces the stochastic time series observed 
in experiments. Stochastic Thermodynamics refers to a 
general framework for a consistent definition of fluctuat¬ 
ing work and heat currents on the level of these fluctu¬ 
ating time series Eli- The common model for molecu¬ 
lar motors are dynamically reversible Markov jump pro¬ 
cesses, which can be thought of as (memoryless) ran¬ 
dom walks on a biochemical network of states mEi]. 


In an accompanying publication m we investigated the 
asymptotic statistics of such systems from the perspec¬ 
tive of the cycle topology of the network of states. In 
particular, we developed an efficient method to calculate 
all cumulants of arbitrary fluctuating currents analyti¬ 
cally. 


Here, we are interested in the first and second order 
fluctuation statistics, i. e. the expressions for macroscopic 
average currents (like the motor’s velocity) and Green- 
Kubo time-correlation integrals (like its diffusion con¬ 
stant). To be concrete, we use the analytic nature of 
our method to analyze the parameter space of different 
stochastic models for the motor protein kinesin mHia, 
which were designed to reflect typical force-spectroscopy 
experiments inHni. Besides illustrating the insights 
that thermodynamic cycles provide into the motor dy¬ 
namics, our results uncover interesting model predictions 
and thus indicate directions for future experimental re¬ 
search: The validity of a nonequilibrium fluctuation dissi¬ 
pation relation at mechanical stalling as well as negative 
differential mobility, commonly referred to as “getting 
more from pushing less” m- 


This work is structured as follows. In Sec. [IT] we briefly 
review the results of Ref. EO]. In contrast to the formal 
exposition there, here we focus on the implementation 
of a universally applicable algorithm for the efficient cal¬ 
culation of averages and correlation integrals of fluctu¬ 


ating currents in Stochastic Thermodynamics. Sec. Ill 


thoroughly discusses how to apply our universal method 
in the concrete biophysical context of a kinesin model. 
In Sec. IV we give a detailed account of kinesin’s chemi- 
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cal (ATP hydrolysis) and mechanical (displacement) cur¬ 
rents as functions of their conjugate chemical and me¬ 
chanical drivings. We conclude in Sec. |V| with a discus¬ 
sion of the main conceptional and biophysical insights. 


II. FLUCTUATING CURRENTS AND THEIR 
STATISTICS 

In this section we introduce our mathematical notation 
and - based on the general results presented in Ref. m 
- provide a concise recipe for calculating the averages 
and asymptotic (co-)variances of two fluctuating currents 
in a dynamically-reversible Markov process on a finite 
state space. Such averages and covariances play a ma¬ 
jor role in Stochastic Thermodynamics [ilinillo]: They 
correspond to physical steady-state currents and time- 
correlation (Green-Kubo) integrals [21]. To be concrete, 
we exemplify topological concepts for both a four-state 
and a six-state Markov process, Fig.[2 In Sec. m we in¬ 
terpret these examples as models for the molecular motor 
kinesin. 


A. Currents for Markovian processes 

Memoryless stochastic processes on a finite state space 
V = {'^ 1 ,'^ 2 , • • • 5 are called Markovian jump pro¬ 
cesses. Henceforth, we consider the time-continuous and 
homogeneous case. A realization, or trajectory (y/ct/^), 
of the process starting at time to = 0 in a state 70 G V 
is a collection of jump times t^ > 0 and visited states 
7 /c G V with /c G N. We interpret it as a time series 
that contains the outcomes of subsequent measurements 
performed on a small system like a molecular motor. 

Transitions from a state G V to a different state 
Vj G V occur at a given constant rate Wj. For thermo¬ 
dynamic consistency IHHl Eol EI], we require dynamical 
reversibility^ i. e. vo'j > 0 4=^ vuj > 0. With this constraint 
we draw the state space as an undirected graph Q with 
the N states as vertices and M admissible transitions as 
edges, cf. Fig. In addition to dynamic reversibility we 
assume that the state space is connected, which ensures 
ergodicity of the process. 

An ensemble of trajectories with initial probability dis¬ 
tribution p( 0 ) = (pi( 0 ),... ,PAr( 0 )) on V evolves accord¬ 
ing to the master equation [ 22 ]: ^p{t) = p{t)W or, in 
components, 

where we use the convention wl = — • Ergodic¬ 

ity of the process implies that there is a unique steady- 
state probability distribution tt satisfying 0 = ttW to 
which all initial conditions will converge eventually. The 
quantities Jj = — iTjwl represent the steady-state 

probability currents between two states Vi and Vj. The 


special condition where all the steady-state currents van¬ 
ish, Jj =0, is called detailed balance or equilibrium. Since 
we are interested in the currents of non-equilibrium sys¬ 
tems, we will not assume detailed balance in the follow¬ 
ing. 

In order to account for general currents, e. g. changes 
in energy, entropy, particle numbers or physical position, 
we introduce jump observables. A jump observable ip 
assigns a weight to the transition from Vi to vj , where 

we require anti-symmetry: cp'j = —(fj. The macroscopic 
average current c((p) associated to a jump observable (p 
is 


hj 


In order to illustrate the concept we provide two ex¬ 
amples. For a pair of states (vi^vj) we define a simple 
but important case of jump observable: the counting ob¬ 
servable (p(ij) counts +1 for a transition from state Vi 
to Vj and —1 for a reverse transition from Vj to Vi. To 
every other transition it associates a weight of 0. In this 
case the macroscopic counting rate from Vi to Vj equals 
the steady-state probability current between these states: 
c{(p(^ij^) = Jj . This expression obviously vanishes if tran¬ 
sitions between Vi and Vj are impossible. In Ref. m we 
emphasized that counting observables form a basis of the 
space of jump observables: every jump observable can be 
expressed as an appropriate linear combination of count¬ 
ing observables. 

Another important example is the dissipation in 
stochastic thermodynamics. It is derived from the jump 

observable a that takes the values a*- = In ^. The 
macroscopic average dissipation 


I w] 

= 2 E ~ 

^ w - 

^,0 * 


is non-negative and vanishes only at equilibrium, i. e. if 
and only if we have detailed balance. 

Kirchhoff’s current law states that the currents in an 
electrical network balance at each vertex vi. The same is 
true for the steady-state probability currents Jj, and the 
stationary Master equation formalizes Kirchhoff’s cur¬ 
rent law as Jj = 0. In Ref. [8] , Schnakenberg dis¬ 
cussed an extended analogy between Markov jump pro¬ 
cesses and Kirchhoff’s laws. The fundamental object 
of Schnakenberg’s theory is a set of H = M — N -\- 1 
fundamental cycles {0}, which describe the topology 
of a Markov jump process. They are obtained from a 
spanning tree of the graph representing the network of 
states, cf. Fig. as well as Sec. El For now, think 
of a fundamental cycle Q as a tuple of consecutive 
states ( 70 , 71 , ••• ,7m(^) = 7o) that form a self-avoiding 
and closed trajectory. Cycles are defined up to cyclic 
permutations. Adding the contributions of a jump ob¬ 
servable cp along the transitions in a fundamental cycle 
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FIG. 1. Two different graphs representing Markov models 
with (a) four states (iV = 4, M = 5) and (c) six {N = 6, 
M = 7) states. The edges marked in green serve as a spanning 
tree T. They connect all vertices of the respective graph. The 
remaining edges, marked blue, are the respective chords Ti. 
Here, we already indicate an orientation for the chords to 
provide a reference for the sign of the currents. Each chord r]i 
gives rise to a fundamental cycle Q . They are shown in panels 
(b) and (d) for the four and six-state models, respectively. 
Regarding the topological cycle structure, both graphs are 
equivalent. In particular, they have the same number B = 
M — iV+l = 2of fundamental cycles. 

Q gives its circulation 

m{£) 

‘fie ■■= E • (2) 

k=0 

Kirchhoff’s voltage law states that the voltage drops 

V = Uj — Ui between vertices of an electric network van¬ 
ish if integrated along any circuit.^ Using the notion of 
circulations, the voltage law reads U = 0. Another result 
obtained in this context is the Schnakenberg decomposi¬ 
tion of the macroscopic dissipation rate [8]: 

B 

= (3) 

e=i 

where q is the steady-state probability current associ¬ 
ated to the cycle Q HO]. The circulations &£ of the dis¬ 
sipation are called the cycle affinities. In the context of 
irreversible thermodynamics [23], the Schnakenberg de¬ 
composition § identifies the cycle affinities as the gen¬ 
eralized forces which are conjugate to the cycle currents 

Cl. 

In an earlier publication [T3j the authors pointed out 
that Schnakenberg’s decomposition is equally applicable 
to other observables ip. As a corollary, one may express 
the average current c{ip) = ci ipi using only the cyclic 
structure of the graph. Fig. shows two graphs with 
four and six states, respectively, but having the same 
cyclic structure. Collecting all B circulations of an ob¬ 
servable in a 5-tuple gives its chord representation 
ip^ := ((^ 1 ,..., (^^) G The detailed mathematical 

background of this representation is discussed in Ref. [ID]. 


Here we are interested not only in the macroscopic 
expectations of currents, but also in their higher order 
statistics. Consequently, the object of study in this work 
are fluctuating currents. The instantaneous current (t) 
derived from a jump observable ip along a trajectory 
( 7 /c, tk) is defined as 

CX3 

k=0 

The time-integrated current 

n{T) 

(4) 

/c=0 

thus accounts for the total change of the observable ip 
along a random trajectory with a random number n{T) 
of jumps up to time T. Hence, this time-integral is a 
random variable with its own statistics. A typical real¬ 
ization of the time-integrated current, and thus its ex¬ 
pectation value, grow linearly in time. Due to ergodicity, 
the time-averaged current ^pj^ := ^ip{T) converges to the 
macroscopic current c{ip) in the long time limit: 

c(^). (5) 

Equivalently, one can average the fluctuating current 
over trajectories in the steady-state ensemble: {jcp{t)) = 
{jip{0)) = c{ip) , where one exploits the fact that the 
steady state is time independent. This ensures that the 
macroscopic current c{ip) is the expectation value, or 
mean, of the fiuctuating current j(^{t). 

Another important statistical measure are the corre¬ 
lations of two currents and j^. A measure for this 
correlation is the Green-Kubo integral: 

POO 

c{<P,'>P) ■= / ((Jv(0) - c{ip)) - citp))) dt. (6) 

Jt=o 

Similar to the case of the average currents, ergodicity al¬ 
lows us to replace the steady-state ensemble average by a 
time average over the argument of the first current j^p. As 
a consequence [21], the correlation integral corresponds 
to a properly scaled covariance of the time-averaged cur¬ 
rents: 

c((/9,'0) = lim T Coy \ipp,f!j^ . (7) 

As such, the macroscopic current and the Green-Kubo 
integral are the first two scaled cumulants of the pair 
of time-averaged currents. Scaled cumu¬ 
lants are defined as derivatives of the scaled cumulant- 
generating function^ which can be obtained using meth¬ 
ods from Large Deviation theory [lOlEIlEllISS]. Higher 
order derivatives represent higher orders of the statistics, 
such as skewness and kurtosis. 

In our accompanying publication Ref. m we prove a 
gauge invariance of the fluctuation statistics and show 
in detail how Schnakenberg’s decomposition is extended 
to all cumulants of arbitrary observables. In the next 
section we give a brief review in form of a concise and 
efficient recipe for the first two orders. 
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B. Determining averages and (co-)variances of 
currents 

The only ingredients needed for the calculation of the 
scaled cumulants are the transition matrix W and the 
jump observables representing the currents of interest. 
The elements rej of the transition matrix as well as the 
jump observables may depend on the (physical) control 
parameters of the Markov processes in an arbitrary way. 
In the following, we present a simple yet efficient algo¬ 
rithm for the calculation of the first two scaled cumulants 
of two jump observables (p and It consists of three sub¬ 
steps: topological, algebraic and physical. The first two 
steps are universal. Only the last step involves the jump 
observables in question. Note that the algorithm does 
neither require the steady-state distribution tt nor the 
scaled-cumulant generating function. We emphasize this 
fact, because in general these quantities are difficult or 
even impossible to obtain analytically, i. e. in the form of 
a fully parameter-dependent, symbolic expression. 


1. Topology: defining fundamental cycles 

The first step in the analysis addresses the topology of 
the graph Q representing a network of states, cf. Fig. 

a. Choose a spanning tree T for the undirected graph, 
i. e. an undirected subgraph spanning all vertices but 
not containing any circuit (green edges in Figs, 
and c). 

b. Provide an orientation to the B = N — M 1 undi¬ 
rected edges r]i eTL left out by the tree T. They are 
the called chords (blue edges in Figs. ^,c). 

c. Identify the fundamental cycles: for every chord r]i G 
H, its terminus and origin are connected by a unique 
directed path through the spanning tree. Adding the 
chord itself as a closure of this path results in the 
fundamental cycle Q (Figs. ^,d). 

2. Algebra: determining the fundamental current cumulants 

The second step of our algorithm involves the deter¬ 
mination of the first two (joint) scaled cumulants of the 
fluctuating currents associated to the chords r]£ eTL. 


c. Calculate the vector c G with entries q = c{r]i) 
and the scaled co-variance matrix C G with 

entries C^m = c{r]£^r]m), as follows: 


C£ 


C£ 


m 


d£ao 

ai 

_ 2{d£ao){dmao)a2 

ai af 

{dmOjl){d£Oj{fj “h {d£Qj\){^dfYi(l{)) 
I 9 


(9a) 


ai 

(9b) 


where the partial derivatives d£ak •= 
the coefficients are evaluated at q = 0. 


and 


Remark: Higher order scaled cumulants are similarly 
accessible. The characteristic equation 0 = Xn{^]Q) 
uniquely defines the entire scaled cumulant-generating 
function \u{q) with A'^(O) = 0. Taking derivatives of 
the characteristic equation yields linear equations for the 
cumulants, i. e. the inner derivatives A'^(O). Note 

that higher order cumulants depend on the coefficients 
cik{(l) with /c > 2, and the symbolic expressions become 
more complex. The first two orders are explicitly given by 
Eqs. (©• The symbolic manipulations that are necessary 
to obtain the higher orders are efficiently implemented in 
modern computer algebra systems. For more details on 
the procedure and a derivation of Eqs. (|^, the reader is 
referred to our accompanying publication m- 


3. Physics: cumulants of jump observables 

The third and final step of the algorithm yields the 
first two scaled cumulants of the fluctuating currents as¬ 
sociated to the jump observables p and fj. 

a. Sum the jump observables p and along the edges of 

the fundamental cycle Q to obtain the circulations 
P£ and f;£. They are the coordinates of the chord 
representations pu , G . 

b. The steady-state average of p^ and of the scaled co- 
variance of p and then read 


a. Write down the characteristic polynomial 
Xw(A; 9 i, 92 , • • • ,te) = det(W« - AI) of the 
matrix with entries 


(w„); = 


w) exp {±qe) 




if (z j) = ±7]i, 

else. 


(8) 


B 

c{(p) = c-(p-H = 'f2‘Pece 
e=i 


B 

c{(p, ^P) = (fu ■ c IpH = Ti Clm tpm 

m,^=l 


(10a) 

(10b) 


b. Identify the coefficients ao{q)^ <^i(q) and a 2 {q) of 
X'h(A;q) = fhe coefficients of 

the constant, the linear and the quadratic term. 


We conclude the section with final remarks on the 
choice of the spanning tree in step 1, which is a pri¬ 
ori arbitrary. Different choices yield different chords - 
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FIG. 2. (a) Kinesin is a motor protein consisting of two iden¬ 
tical entangled subunits. Cargo is bound at the tail end. The 
active sites on kinesin’s head end bind to the microtubule and 
act as kinesin’s “feet”, enabling the molecule to perform di¬ 
rected steps. Kinesin’s stepping mechanism is the result of 
subsequent changes in how strong the active sites bind to the 
microtubule. The trailing (left) and leading (right) sites are 
represented by colored ellipses. ATP-laden (red) and empty 
(gray) sites bind strongly, whereas an ADP-laden (blue) site 
binds only weakly. The succession of chemical compositions 
shown in panel (b) is called the forward cycle T. Starting from 
the upper left state, the forward cycle involves (in counter¬ 
clockwise direction): (i) Binding of ATP to the (empty) lead¬ 
ing site (ii) a mechanical step (brown edge), i. e. the exchange 
of the leading and trailing site, (iii) release of ADP from the 
(new) leading site, (iv) hydrolysis (purple edge) of ATP into 
ADP at the trailing site. 


in alternating succession, thereby allowing the motor to 
perform mechanical steps of length L = 8nm [301131]. 
Due to the polarity of the microtubule, this motion has 
a preferred “forward” direction. 

The energy necessary for this active directed transport 
is provided by the hydrolysis of adenosine triphosphate 
(ATP) into adenosine diphosphate (ADP) and inorganic 
phosphate (P). Unlike macroscopic motors, small molec¬ 
ular machines operate at low Reynolds numbers and in¬ 
ertia plays no role: chemical energy is not converted into 
mechanical energy by a transfer of momentum. Instead, 
kinesin’s mechanical displacement is the result of a com¬ 
plex interplay of the strength of the microtubule bind¬ 
ing at the active sites, which depends on their chemical 
composition. ATP-laden and empty sites bind strongly, 
while ADP-laden sites bind weaker [igin]. Under physi¬ 
ological conditions the mechano-chemical interaction can 
be described by the “forward cycle” depicted in Fig. 
HU. Models that only treat the forward cycle feature 
tight coupling between the hydrolysis reaction and the 
stepping: each hydrolysis of an ATP molecule gives rise 
to exactly one motor step m- 


B. Experiments and models 


and thus different expressions for the fundamental cur¬ 
rent vector c and the fundamental covariance matrix C. 
Expressions and ( [1Q| ) are universal. In order to cal¬ 
culate the cumulants of any jump observable, no equa¬ 
tions need to be solved. Via Eq. any combinatorial 
complexity is hidden in (the derivatives of) the coeffi¬ 
cients Ok of the characteristic polynomial. The latter are 
calculated in a straightforward way either manually or 
by using a computer-algebra system. However, the final 
expressions ([Tq|) have fewer terms if some of the circu¬ 


lations (p£ or 'ipi along fundamental cycles vanish. It is 
thus worthwhile to take a careful look at the particular 
set of jump observables cp and under consideration and 
choose a spanning tree that is optimal in that regard. 


III. KINESIN 

A. Kinesin as a molecular motor 

Kinesin is a molecular motor which facilitates trans¬ 
port in eukaryotic cells. It moves along intracellular fila¬ 
ments called microtubules and plays a major role in many 
biological processes, including mitosis, meiosis and trans¬ 
port of cellular cargo. The most well studied variety of 
kinesin - both experimentally (see e. g. 0 HSHni and 
References therein) and theoretically [TTl [121 (SMUj “ is 
a protein dimer consisting of two identical subunits. Eig- 
ure shows a sketch of kinesin binding its intracellular 
cargo at its tail end. Kinesin’s head end consists of two 
active sites which bind and unbind to the microtubule 


An important biophysical question regards the force 
that kinesin generates for different concentrations of the 
chemicals ATP, ADP and P involved in the hydrolysis 
reaction. Typically, experiments measure this force by 
linking kinesin to a di-electric colloidal bead which resides 
in an optical trap [HHISlEo]. Involved experimental set¬ 
ups allow the precise control of the pulling force F that 
the optical trap exerts on the motor against its typical 
direction of motion. The independent driving parame¬ 
ters are the non-dimensionalized force / := {LF)/{kBT) 
and the non-dimensionalized chemical potential differ¬ 
ence A/i = log (iFeq[ATP]/[ADP][P]), where /cp denotes 
Boltzmann’s constant, T the temperature, Keq the equi¬ 
librium constant of the hydrolysis reaction and [X] the 
concentration of chemical species X. In the remainder of 
this work, all physical quantities are expressed in units 
based on the length scale L, time scale 1 s, and energy 
scale k^T. 

Many experiments probe the stalling force /stall 
which is defined as the value of the force needed to bring 
the motor to a halt for a given chemical potential dif¬ 
ference A/i. Under physiological chemical conditions, 
kinesin hydrolyzes ATP even at stalling forces [mnn. 
The exact details of the kinesin stepping mechanism un¬ 
der high mechanical loads remain unknown and several 
models exist, cf. Refs. HU Ha HZ] and the references dis¬ 
cussed in these publications. While these models differ in 
their details, they all feature more than only the tightly 
coupled forward cycle. 

A prominent example of a thermodynamically consis¬ 
tent model was introduced in Ref. HU. There, the key 
idea is to extend the forward cycle shown in Eig. It by 
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FIG. 3. (a) The six-state kinesin model from Ref. m extends 
the forward cycle from Fig. 1^ by two states. For this model, 
We use the same spanning tree and chords as in the example 
shown in Fig. [^. Then, the fundamental cycles and 

C ,2 = 'D are the forward and dissipative cycle (b) and (d), 
respectively. The backward cycle (c) is the linear combination 

R = C2 - Cl, c/. Ref. [To]. 


the chemical compositions obtained from exchanging the 
trailing with the leading active sites, cf. Fig. [^. In ad¬ 
dition to the forward cycle T^ the extended network fea¬ 
tures six states and has two additional cycles, Figs. |2j3-d: 
the backward cycle B represents backward motion under 
hydrolysis, while the dissipative slip cycle V represents 
the futile hydrolysis of two ATP molecules without any 
stepping [27]. In such a multiple-cycle model, hydroly¬ 
sis and mechanical displacement are not tightly coupled 
anymore. In Sec. |IV A| we will address the question of 
quasi-tight coupling^ i. e. situations where the ratio of the 
average number of chemical and mechanical events pre¬ 
dicted by the model is close to unity. 


C. Network theory for the kinesin model 

In order to study quasi-tight coupling, energy conver¬ 
sion and the predicted response to changes in the driving 
parameters, we apply the algorithm presented in Sec. EH 
to the six-state model for kinesin. Fig. For the first 
step of the algorithm, we choose the spanning tree and 
its chords in the same way as in Fig.[^,d. Consequently, 
the fundamental cycles C,i = T and (2 = correspond 
to the forward and dissipative cycles, respectively. 

The second step of the algorithm requires the deter¬ 
mination of the fundamental current vector c and the 
fundamental co-variance matrix C. With the enumera¬ 
tion of the vertices as in Fig. the matrix W'^(gi,g 2 ) 
reads: 


/^i 

wl 

0 

0 

0 

wl\ 

Wi 

wl 

wl 

0 

wl 

0 

0 

wl 

wl 

wl 

0 

0 

0 

0 

wl 

wf 

wie^‘^ 

0 

0 


0 


wl 

wl 

\wf 

0 

0 

0 

wl 

wlJ 


It is straightforward to write down its characteristic poly- 
nomial x«(A; 91 ,^ 2 ) =: Yfk=o 5 ' 2 )A*, and to extract 


the coefficients ao{q) = det W^(gi, ^ 2 ), cii and a 2 . Dif¬ 
ferentiating with respect to qi and g 2 , and evaluating at 
Qi = Q 2 = ^ yields the expressions d^ak appearing in 

Eqs. 0. 

The third step requires the circulations of the jump 
observables of interest. For the present discussion, we 
consider the displacement d = ^( 2 , 5 ) and the hydrolysis 
count h = ^( 6 , 1 )+^( 3 , 4 ), which indicate a transition along 
the brown and purple edges in Fig.[^ respectively. Their 
matrix representations read 

d] = (11a) 

h] = (Hb) 

where denotes the Kronecker delta, which yields 

one if m = n and zero otherwise. The circulations of 
d and h simply count the number of the brown and pur¬ 
ple edges in the fundamental cycles Ci = T and C 2 = 
cf. Fig. |3]3 ,d. The chord representations of d and h thus 
are 

du = {di, d 2 ) = (1,0), (12a) 

hH = {k,h2) = {l,2). (12b) 

Note that the choice of the chords is optimal for the 
calculation of the present variables because one of the 
entries of dy^ vanishes (^2 = 0), while this cannot be 
achieved for h^. After all, all cycles contain at least one 
hydrolysis event. 

The corresponding macroscopic currents, i. e. the ve¬ 
locity c{d) and the hydrolysis rate c{h) are obtained from 
Eq. ( p^ as: 

c{d) = Cl, (13a) 

c{h) = Cl + 2 c 2 . (13b) 

Their scaled (co-)variances amount to 
c{d,d) = Cii, 

c(/i, h) = + 4(7i2 + 4(722, 

c(/i, d) = (7ii + 2(7 i2. 

In addition to displacement d and hydrolysis count 
/i, we are interested in the jump observable crj = 

ln{wj/wl) corresponding to the dissipation. As dis¬ 
cussed in Sec. EH its circulations are the cycle affini¬ 
ties. The Hill-Schnakenberg conditions are necessary 
for the consistency of a Markov jump process with the 
thermodynamic notion of local equilibrium |3l|9]. They 
state that the affinity of a cycle must express the (non- 
dimensionalized) differences in the potentials of the reser¬ 
voirs, cf. Refs. [aiHiiin]. Upon completing the forward 
cycle T = (i, an amount A/i of chemical energy is used 
by the system to perform a (dimensionless) amount —/ of 
work against the pulling force. Similarly, a completion of 
the dissipative cycle V = (^2 uses 2A/i of chemical energy. 
Consequently, the chord representation of the dissipation 
reads 


cr-H = {di, (72) = (-/ + A/i, 2A/i). 
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(a) . (b) (c) (d) 

d \ f * t0t 

^ 4P ^ dP^— ^ 


FIG. 4. A model with four states, which describes the same 
physics as the six-state model shown in Fig.[^ In this simpler 
model we combined the transitions for the ATP hydrolysis on 
one active site with the ADP release on the other one into a 
single transition. Details of the model construction are given 
in App.[^ 


The Schnakenberg decomposition thus lets us express the 
average steady-state dissipation by physical parameters 
and currents through fundamental chords 

c(o') = (-/ + A/i)ci + (2A/i)c2. (14) 


D. The cycle perspective 


In the previous section, we expressed observable quan¬ 
tities only by means of their circulations around funda¬ 
mental cycles and the fundamental first and second chord 
cumulants. Nowhere in these expressions do the num¬ 
ber of states or the choice of a spanning tree appear ex¬ 
plicitly. Hence, the same expressions are reproduced by 
any model with the same cycle topology - as long as the 
physics along the cycles, i. e. the circulations of antisym¬ 
metric jump observables, are the same. As exemplified 
in Sec IIA in Fig. [^,b, one can formulate a model on 
four states with the same cycle topology as the six-state 
model described in Fig.|^ Allowing only for single edges 
between states, such a four-state model is the minimal 
model featuring two independent cycles. An interesting 
question is how this model (and other reduced models) 
compare to more complicated ones. 

In Ref. [13] we used the idea of preserving the cy¬ 
cle affinities and circulations of jump observables along 
cycles (together with locality constraints) to develop a 
coarse-graining algorithm for stochastic models. Its ap¬ 
plication to the six-state kinesin model produced various 
topologically equivalent models, which all preserved the 
fluctuation statistics of the observables of interest almost 
perfectly. A disadvantage of this coarse-graining algo¬ 
rithm is that the individual transitions in the network of 
states lose their original interpretation. 

In contrast, Fig. shows a four-state model with a 
clear interpretation of the transitions. This has the ad¬ 
vantage that the parameterization of the transition rates 
is found by the same physical arguments as the ones used 
in Ref. m for the six-state model. Details of the con¬ 
struction of this model are given in Appendix We will 
see in the next section that all the predictions of the 
six-state model are also found in topologically equiva¬ 
lent four-state models. This observation underlines the 


virtue of viewing models in ST from the perspective of 
cycles - an idea that was pioneered by Hill 13 HH and 
Schnakenberg |8] and has regained considerable attention 
recently [IIl[I3[271|32ll33]. 


IV. RESULTS 

One of the main messages of this work is that the algo¬ 
rithm presented in Sec. |HB| allows us to probe parametric 
models used in Stochastic Thermodynamics in a system¬ 
atic way. In order to demonstrate the efficiency of our 
approach, we report on various non-trivial predictions of 
the six-state kinesin model. At the end of this section we 
will compare these results with other models. 

Throughout this section, we choose the parameter 
range similar to Ref. [27| and vary —30 < /, A/i < 30. 
Then, in physical units the pulling force F varies between 
about —15 and 15 pN. Following Ref. m, the chemical 
potential difference is adjusted by changing the ATP con¬ 
centration while fixing the other chemical concentrations 
at physiological values (see also Appendix]^. The physi¬ 
ologically relevant region for the chemical driving param¬ 
eter is limited to about 20 < A/i < 30. Negative values of 
the chemical potential correspond to extremely low ATP 
concentrations. In particular, the “homeopathic limit” is 
reached at about A/i < —14: At that point there is less 
than one ATP molecule in an experiment containing one 
liter of solution. 

The reason we still present our results for the entire 
parameter range —30 < /, A/i < 30 is two-fold: firstly, it 
enables a direct comparison with previous work [niEzi. 
Secondly, it demonstrates the effectiveness of our algo¬ 
rithm in predicting results that vary over many orders of 
magnitude. Still, we emphasize that non-trivial results 
are encountered exactly in the experimentally accessible 
region, where we consider the model as valid. 


Velocity, hydrolysis rate, and their quasi tight 
coupling 


Fig. reproduces a central result of Ref. m con¬ 
cerning the operations modes of kinesin. These modes 
are defined by the signs of the average currents, i. e. of 
the velocity c{d) and the hydrolysis rate c{h). How¬ 
ever, the resulting phase diagram contains no informa¬ 
tion regarding their magnitude. Based on the expres¬ 
sions (13) we provide a detailed account on their numer¬ 
ical values in Fig. If- Note that these currents vary over 
about 20 orders of magnitude. This underlines the im¬ 
portance of having analytical expressions to generate the 
plots. A brute-force numerical approach will either be 
prohibitively expensive in terms of computer resources, 
or it will suffer from severe inaccuracies when dealing 
with this vast range of numerical values. 

The analytical expressions for the currents also reveal 
an interesting relation between the average currents. In 
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(a) operation modes (b) velocity c(d) hydrolysis rate c(h) (c) current ratio c{h)/c{d) 



FIG. 5. (a) Operation modes as identified in Ref. [27]. (b) Decadic logarithm of the absolute value of the average velocity 
c{d) (left) and hydrolysis rate c(h) (right). The white lines indicate where the currents vanish, i. e. the solid lines displayed in 
panel (a). The contour lines show that the macroscopic currents are proportional away from these lines, (c) Plotting the ratio 
c(h)/c[d) makes this proportionality visible directly. The proportionality constant has an absolute value very close to unity, 
indicating quasi-tight coupling for most parameter values (see discussion in the main text). Note thate the values of the ratio 
are cropped at absolute values of 2, most prominent in the dark regions surrounding the singularities of the ratio. 


Fig. we plot the ratio c{h)/c{d) of the hydrolysis rate 
and the velocity. Again, access to the analytical expres¬ 
sions for the currents is crucial to determine the ratio. 
After all, both its numerator and denominator are of 
the order 10“^^ in the lower left corner of the param¬ 
eter space. 

The most prominent feature of Fig. HI is that away 
from the zero-current lines, the ratio of average hydroly¬ 
sis rate and velocity takes values very close to ±1. Con¬ 
sequently, on average the completion of a cycle yields 
one mechanical step and one chemical event. We say 
that chemical and mechanical currents are quasi-tightly 
coupled. Experimentally, it was found that kinesin hy¬ 
drolyzes one ATP molecule for each mechanical step El- 
According to the model considered here, quasi-tight cou¬ 
pling is a generic feature that holds more generally: even 
in the region where kinesin moves backward while con¬ 
suming ATP, the absolute values of the currents are 
locked to a ratio of one. 

Knowing the absolute values of the currents rather 
than only their signs also allows us to treat kinesin’s 
thermodynamic cycles in more detail. In Ref. this 
discussion was based on the signs of Hill’s (excess) cy¬ 
cle fluxes [7]- With the current ratio we interpret the 
regions shown in the phase diagram Fig. in terms of 
dominant cycles - at least away from their boundaries: in 
the upper left and lower right regions the forward cycle 
dominates such that average hydrolysis and velocity 
are directly proportional, c{h) ^ c{d). The difference 
between those regions is the angular direction: counter¬ 
clockwise completion leads to a forward movement ac¬ 
companied by ATP hydrolysis, whereas clockwise com¬ 
pletion yields backward stepping and ATP synthesis. In 
contrast, in the upper right and lower left regions hydrol¬ 
ysis and velocity are anti-proportional, c{h) c{d): a 
counter-clockwise (backward, hydrolysis) or a clockwise 
(forward, synthesis) completion of the backward cycle B 


dominates the average dynamics, respectively. This re¬ 
sult, which is based on the values of physiological cur¬ 
rents, thus complements and extends the discussion pre¬ 
sented in Ref. m- 

B. Efficiency of energy conversion 

Under physiological conditions, kinesin uses the chem¬ 
ical energy released by the ATP hydrolysis to perform 
mechanical work. Energy efficiency is one of the most 
important questions for molecular machines involved in 
cellular energy conversion [34f[36] . just as it is for macro¬ 
scopic machines. A framework for a quantitative analysis 
is based on the notion of conjugate currents and forces 
from irreversible thermodynamics [23] . Generally, a com¬ 
plete set of conjugate currents c{(pi) and forces Ei yields 
the average dissipation as the bilinear form 

i i 

where = (piEi denotes the distinct contributions to 
the entropy production. 

Using equations © and ( pAl we find that 

c(cr) = {-f)c{d) + {Aia)c{h) 

—. c((J]xiech) “b c((Jchem)* 

We see that in the kinesin model, velocity c{d) and hy¬ 
drolysis rate c{h) are conjugate to the negative pulling 
force — / and the chemical potential A/i, respectively. 
Eor a system with two independent contributions to the 
entropy production, a = cri -|- <72, one may define the 
efficiency of energy conversion in general terms [35]. To 
that end note that c(cr) is always positive. This, however, 
does not imply that both contributions c(cr^) are positive. 
Indeed, systems act as energy converters only if one of 
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the contributions, say di, is negative. Then, a (positive) 
average power output Wout •= ~c(cri) is sustained by a 
(positive) average power input Win = <^((72). Note that 
c{(j 2 ) is positive and larger in magnitude than c(cri), be¬ 
cause c{a) = c(cri) + c(cr2) > 0 must always hold. Hence, 
the efficiency of energy conversion is defined as 


It is always positive and smaller than unity. 

In the framework of Stochastic Thermodynamics, this 
efficiency has been studied under various aspects (c/. 
e. g. Refs. [34ll37] b In Fig. we give the efficiency of 
energy conversion 7 ) for the kinesin model. The regions 
A-D correspond to different types of energy conversion 
where the system either acts as a motor (A,C) or a chem¬ 
ical factory (B,D). Outside of these regions both contri¬ 
butions to the entropy production are positive and no 
energy conversion takes place. 

We note the following prediction: for any fixed value 
of A/i in the physiological range, i. e. for 20 < A/i < 30, 
the value of the force at maximum efficiency is around 
/ ^ 10.5. This suggest that kinesin might be optimized 
to encounter (elastic) forces of around 5 pN, independent 
of the ATP concentration. It will be interesting to ex¬ 
plore the implications of this result for the collaborative 
behavior of multiple kinesin molecules involved in the 
viscous transport of organelles. 

Finding the parameters of a system that extremize 
thermodynamic quantities is a generic problem. Re¬ 
cently, many authors have discussed the notion of effi¬ 
ciency at maximum power (see e. g. Refs. [34l [38] and 
therein). Having fully parameter-dependent symbolic ex¬ 
pressions for the various contributions to the entropy pro¬ 
duction establishes a general (analytic) approach to this 
optimization problem. 


(^2 


< 1 . 


(16) 
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0.2 
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FIG. 6. Efficiency of energy conversion in the kinesin model. 
The four regions A-D correspond to four different ways of 
energy conversion. In the regions outside the solid curves, 
no conversion between mechanical and chemical energy takes 
place. In regions A and C kinesin acts as a motor converting 
chemical into mechanical energy against the external force. 
In regions B and D kinesin resembles a chemical factory that 
uses mechanical energy to produce ATP and ADP, respec¬ 
tively, against the chemical potential provided by the solu¬ 
tion. The sketches in the upper right and lower left illustrate 
the combination of thermodynamic forces acting on the mo¬ 
tor in the respective quadrant. In the upper right, kinesin is 
pulled backward (z. e. against its standard direction of mo¬ 
tion) in an ATP rich environment. In the lower left, kinesin 
is pushed forward in an ADP rich environment. Energy con¬ 
version only occurs in the regions where both mechanical and 
chemical currents have the same sign and the forward cycle 
dominates, c/. Eig.[^ 



C. Diffusion constant and randomness parameter 


So far, we have only investigated average currents, 
which are also available if the steady-state distribution 
TT is known, cf. Eq. 0. Higher order statistics of fluc¬ 
tuating currents cannot be expressed by means of the 
stationary distribution only, although a perturbation ex¬ 
pansion exist [39|. The method presented here provides 
direct ac cess to the (co-)variance of fluctuating currents 
via Eq. (10b), without the knowledge of the stationary 


distribution. 

For motor proteins we are mostly interested in the sec¬ 
ond scaled cumulant of the time-averaged displacement. 
It quantifies the (linear) scaling of the (fluctuating part 
of) the mean-square displacement, and thus defines (up 
to a factor of two) the non-equilibrium diffusion constant 
D = D(/,A/i): 


d) = ^ {{d{Tf) - {d{T)f^ 


=: 2D. 


(a) Diffusion constant 

D — - c{d, d) 



(b) Inverse randomness parameter 
^ c(d) 
c{d, d) 


-30-20-10 0 10 20 30 

/ 
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FIG. 7. (a) The diffusion constant D = |c(d, d) on a 

decadic logarithmic scale, (b) The inverse randomness pa¬ 
rameter = c(d)/c(d, d) = c{d)/(2D) compares the veloc¬ 
ity and the diffusion constant. Away from the stalling line, 
it obtains an absolute value close to unity. Solid lines show 
where |r| = 1 holds exactly. 
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In Fig. we show the diffusion constant in the six- 
state kinesin model. Like the average velocity, its values 
span a range of about 20 orders of magnitude. Under 
physiological conditions (/ = 0, A/i ^ 25) the diffusion 
constant is about ten orders of magnitude larger than at 
equilibrium. 

A direct measurement of the parameter dependence of 
D is difficult. An observable that is more easily accessi¬ 
ble in experiments is the so-called randomness parameter 
(sometimes called Fano factor) [151 HZl HOI SI] 


r = lim 

T^oo 


{d{Tf) - {d{T)f 
{d{T)) 


c{d, d) 
c{d) 


(17) 


It is a dimensionless measure of the temporal irregularity 
of the mechanical displacement. While r = 0 indicates 
a deterministic motion without any ff net nations, a value 
of |r| = 1 amounts to a Poisson motor |40]. In Fig. 
we plot its inverse, which is a smooth function. We 
see that the six-state model predicts Poissonian behav¬ 
ior |r| ~ 1 in a large area away from the stalling lines. 
This is in agreement with recent experimental results and 
theoretical predictions from an alternative model m- 
Our method to calculate the second scaled cumulant 
and thus the diffusion constant avoids all of the com¬ 
binatorial complexity of previous approaches mi HD. 
Ref. [43] treats drift velocity and diffusion in Markovian 
models formulated for a periodic lattice in arbitrary di¬ 
mensions. In the present work the topology of physical 
space is independent from the structure of the graph, 
which represents the topology of the model: if a system 
like a molecular motor moves in more than one spatial 
dimension, one defines a distinct jump observable di for 
each of these dimensions i. Up to a factor of two, the 
scaled covariance matrix c{di^dj) then equals the diffu¬ 
sion tensor. 


D. Response theory 
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FIG. 8. Normalized mechanical response Tmech- On the 
dashed green curves, a Green-Kubo FDR ( [T^ holds. One 
of these curves coincides with the stalling line / = /stall (A/i) 
where the average velocity vanishes (white line). The dashed 
blue line indicates a vanishing transport coefficient. To its 
right lies a region of negative differential response. 


cumulants c{(pi^ ipj). The fiuctuation relation for the en¬ 
tropy production [iiio] ensures the validity of the fol¬ 
lowing equilibrium FDR [21] [24] : 

With analytical expressions for the average currents 
c{(pi) one can calculate their derivatives Rij. Because 
the correlation integrals c((/2i, (pj) are known, our method 
enables us to probe the non-equilibrium response proper¬ 
ties predicted by models of Stochastic Thermodynamics. 
As an example, we discuss kinesin’s mechanical response 
using the normalized response coefficient 


Eq. (15) states that the average velocity c{d) and hy¬ 


drolysis rate c{h) are conjugate to the mechanical and 
chemical driving forces — / and A//. Response theory 
studies the dependence of averaged currents J = (c(pi))^ 
to the conjugate fields E = (E^)-. For B independent 
conjugate current-field pairs, {c{pi), Ei)i, the response 
matrix R(-®) is 8i B X B matrix with entries 


^iAE) := 


dc{tpi 


dEi 


(18) 


Fluctuation dissipation relations (FDR) relate the re¬ 
sponse of average currents to their fiuctuation statistics. 
In particular, the Einstein relation relates the mobil¬ 
ity of a particle (or its inverse, the friction coefficient) 
to its diffusion constant. So called Green-Kubo rela¬ 
tions [44] express equilibrium transport coefficients by 
time-correlation integrals, Eq. Here, these time- 

correlation integrals are obtained as second-order scaled 


Tmech (/: A/i) : — 


2 dc{d) 
c{d, d) d{-f) 


1 dc{d) 


( 20 ) 


The equilibrium EDR ( [l^ implies that Tmech(0,0) = 
1. As the transition matrix depends smoothly on the 
driving parameters (—/, A/i) [11], we expect that there 
will be a one-dimensional curve in parameter space where 

Tmech (/: A/i) = 1. 

Eigure depicts Tmech- As expected, we see that 
Tmech(/: A/i) = 1 holds aloug two lincs originating from 
the origin, such that a nonequilibrium EDR holds for 
these parameter values. Remarkably, one of these lines 
coincides with the stalling line / = /stall(A/i), i. e. for 
parameters where the average velocity vanishes. 

Another non-trivial feature of Eigure is the region 
where the normalized mechanical response is negative. 
Since the diffusion constant D is positive, Tmech < 0 
implies that the derivative dc{d)/d{—f) of the mechani¬ 
cal current with respect to its conjugate force is negative. 
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FIG. 9. Comparison of our four-state model with the six-state 
model of Ref. HU. We show the relative errors X^jX^ — 
1 of the corresponding quantities X^ and Xq calculated in 
the four- and six-state models, respectively. Throughout the 
parameter range considered, they are almost everywhere well 
below 15%. Note that for the average hydrolysis rate the 
relative error diverges close to the line c(h) = 0 where the 
hydrolysis rate in the six-state vanishes. Remarkably, this is 
not the case for the average velocity, where the stalling lines 
in both models agree exactly. 


This phenomenon is known as negative differential mobil¬ 
ity [45], or more generally, negative differential response 
(NDR) [T8l[46]. The kinesin model predicts negative dif¬ 
ferential mobility for large enough pulling forces beyond 
stalling, i. e. in situations where the motor walks back¬ 
wards. Then, by pulling more one gets less, i. e. the ve¬ 
locity in pulling direction becomes smaller. This feature 
might already be visible in the experimental data found 
in Refs. [SlIIT]. Although we do not expect to see NDR 
for arbitrarily high pulling forces in real experiments, ex¬ 
plicitly looking for it in the region for small superstalling 
forces seems worthwhile. 


E. Model comparison 

Direct access to the non-trivial features of physical cur¬ 
rents allows us to compare different models in detail, 
both qualitatively and quantitatively. We start by quan¬ 
titatively comparing the results for the six-state model 
(Fig.i with the simpler four-state model (Fig.[^. Recall 
that the latter is constructed following the same physical 
arguments as the former (c/. Appendix]^ for the details). 
The results plotted in Figs. [5]-[^ are all derived from c((i), 
c{h) and c{d^ d) = 2D. In Fig. we plot the relative 
deviations of these quantities between the four-state and 
the six-state model. Throughout most of the parameter 
space, they are only a few percent. This is remarkable, 
because the observables themselves vary over many or¬ 
ders of magnitude. Note that at the boundaries between 
different operation modes (Fig. |^), the first cumulants 
vanish. Hence, we have a divergence in the relative errors 
unless this happens exactly at the same parameter values 
in both models. 

For the hydrolysis rate c{h) such a divergence is vis¬ 
ible in Fig. around (/, A/r) (16,14). In principle, 
this divergence is present wherever c{h) vanishes in the 
six-state model. In practice, however, the curves of zero 


average hydrolysis rate c{h) = 0 agree almost perfectly, 
such that the region where the divergence has an effect 
is extremely small. For most parameters it is hidden 
due to the finite plotting resolutions, and thus not visi¬ 
ble in Fig. [^. In contrast, the prediction of the stalling 
forces /stall (Am) agrees exactly between the two models: 
Fig.|^ does not exhibit any singularities. In Ref. [T3| we 
introduced a coarse-graining procedure which preserves 
the cycle topology of a model. By construction, the first 
cumulants of all currents agree between the original and 
the coarse-grained models. Moreover, the relative error 
in the diffusion constant is comparable in magnitude to 
what we see in Fig. These quantitative results em¬ 
phasize the value of the cycle perspective introduced in 
Sec |IIID[ In order to construct thermodynamically con¬ 
sistent models, one should think of the physics of cycles 
rather than focusing only on individual transitions. 

Finally we compare the six-state model to a general 
model for molecular motors presented in Ref. @3, which 
was studied in detail in Ref. m- Unlike the six-state 
model studied fo far, that model features only two states, 
which correspond to a strongly and a weakly bound con¬ 
figuration. Multiple transition between these two con¬ 
figurations are possible and represent either an active 
(z. e. accompanied by a chemical event) or passive dis¬ 
placement along the microtubule. The cycles of the two 
models are different both in their topology and their 
interpretation. In particular, the two-state model of 
Refs. [HJ [47| has no reference to the “hand-over-hand” 
stepping mechanism of the forward cycle of Ref. m, de¬ 
picted in Fig. IT- Moreover, the two-state model was 
fitted to the experiments of Refs. [HHS], whereas the 
six-state model used the experimental data from Ref. m 

Due to the simple structure of the two-state model, an 
analytical parameter-dependent expression of the scaled 
cumulant generating function was found in Ref. [12]. 
Consequently, analytical expression for the scaled cumu¬ 
lants are known and can be compared to the results ob¬ 
tained for the six-state model of Ref. m, which we used 
so far. Due to the different nature of the models, we do 
not expect their predictions to agree quantitatively. In 
particular this is the case for parameter values that are 
far away from values that are realized in the actual exper¬ 
iments, i. e. for small (or even negative) chemical poten¬ 
tials A/i, or negative values of the pulling force. However, 
for experimentally accessible parameters, it makes sense 
to look for qualitative agreements in the features of the 
two different kinesin models. 


Fi g. [To] shows the normalized mechanical response 
Eq. (20) in both models for sensible chemical potential 
differences (5 < A/i < 30) and positive pulling forces. 
As expected, the models do not agree quantitatively. In 
particular, the stalling lines are at different positions. 
However, they show the same qualitative features: the 
validity of an Einstein EDR at stalling, and the emer¬ 
gence of negative differential mobility just above stalling. 
Together with the experimental hints from Ref. [16], we 
consider this agreement as evidence that negative differ- 
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FIG. 10. Normalized mechanical response Tmech for the six- 
state model of Ref. m studied in the present work (left), and 
the model from Ref. m (right). For experimentally sensible 
parameters both models predict the same qualitative behav¬ 
ior: The validity of a non-equilibrium Einstein FDR (green 
curves) at stalling (white curves) as well as negative differ¬ 
ential response (regions to the right of blue curves) for su¬ 
perstalling forces. Other features (like the overall structure, 
magnitude of the response) show the same qualitative behav¬ 
ior of the two models. 


ential mobility is a generic feature of kinesin — a predic¬ 
tion that should be studied by future experiments. 


V. CONCLUSION 

In the present work, we gave an explicit procedure for 
the analytical, i. e. fully parameter-dependent, calcula¬ 
tion of the statistics of fluctuating currents in Stochastic 
Thermodynamics. The algorithm applies to any finite 
Markov model. We focused on its efficiency in explor¬ 
ing the parameter space of models for the motor protein 
kinesin, while the mathematical background of the algo¬ 
rithm was the subject of an accompanying paper m- 
In the following we summarize conceptual and physical 
insights. 

From a conceptual point of view, we find the following 
points particularly noteworthy: 


Eq. (20) and approach (thermodynamic) optimiza¬ 


tion problems {e. g. efficiency at maximum power, 
cf. Sec.lrV^. 


From a physical perspective, our method allows the 
systematic comparison of the predictions made by var¬ 
ious kinesin models, cf. Sec. |III[ In particular, we gave 
a detailed account on (quasi-)tight coupling, efficiency 
of energy conversion, diffusion and mechanical response 
for a well-known kinesin model [U in Secs. [TVAH iV^ 


Moreover, in Sec. |IVE| we compared these predictions 
with other models. Regarding the modeling of molecular 
motors, we emphasize the following: 


• Current statistics correspond to experimentally ob¬ 
servable quantities, like the average motor velocity 
or its nonequilibrium diffusion constant. Our sys¬ 
tematic approach thus extends and unifies previous 
approaches for calculating these quantities [4TH43] . 


• Thinking of stochastic models in terms of its phys¬ 
ical cycles is useful. It allows model reduction, cf. 
Sec. IV E and Ref. [13]. 


• Independent models predict two interesting 
nonequilibrium response properties of kinesin: 

(i) the validity of a non-equilibrium fluctuation- 
dissipation relation at mechanical stalling, and 

(ii) negative differential mobility for superstalling 
forces. Both of these predictions lie in realistic 
parameter regions and can be tested in future 
experiments. 


Modeling the dynamics of molecular motors as random 
transitions on a biochemical network of states is only one 
of many appliations of finite Markovian jump processes. 
The methods established in the present paper apply to 
any other dynamically reversible model and are easily 
extended to systems with multiple transitions between 
states. Thus, they provide a powerful framework to fully 
explore the physical predictions of any model described 
by Stochastic Thermodynamics. 
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The authors thank Nigel Goldenfeld, Matteo Polettini 
and Massimiliano Esposito for insightful discussions. JV 
acknowledges a research grant of the “Center for Earth 
System Research and Sustainability” while the final ver¬ 
sion of this manuscript was drafted. 


Appendix A: Transition rates of the four-state 
model for kinesin 

The parametrization of the kinesin model on four 
states (Fig. [3 follows the steps in Ref. [TT] for the six- 
state model ^ig.|^. Transition rates 
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Mechanical transition 
Chemical transitions 

(forward cycle) 

Chemical transitions 

(backward cycle) 
Mechanical load 


kI = 3x 10^ 

K.I = 100 

4 = 2.52 X 10® 

= (#)Ci=6.4xlO-“ 
= 2.52 X 10® 


kI = 0.24 

4 = 2.0 

^ 1^3 

4 = 4 = 2.0 
K2 = 4 = 49.3 
Xi = xi = xi = xi = 0.25 


K2 

x! = X3 = X 2 = Xi=0.15 


TABLE 1. Numerical values of the parameters of the four-state model for kinesin. All first-order reaction rates k are given in 
units of s“^ or, if attachment of n chemicals is involved, s“^pM“^. Values correspond to the experimental data in Ref. m as 
stated in Ref. m- The equilibrium constant of the ATP hydrolysis reaction is = 4.9 x 10^^ pM. The parameter 0 = 0.65 
enters the mechanical factor of the transition rates. 


are obtained as first-order rate constants 4 , which are 
multiplied by concentration and force-dependent factors. 
In accordance with first-order rate kinetics, the chemical 
factor reads 

if compound X is attached. 

1 else. 

For chemical concentrations which are not too high, 
the non-dimensional chemical potential difference is given 

by A/u = In (i^eq [adphp] ) ’ where K^q = 4.9 x 
is the chemical equilibrium constant for the ATP hydrol¬ 
ysis reaction happening at kinesin’s active sites. Like in 
Ref. m we fix [ADP] = [P] = 1 pM at physiological 
values and consequently vary the concentration of ATP 
as 

[ATP] = — 

The force dependent factors depend on the non- 
dimensionalized pulling force / = LF/{kBT). They have 
a different form for mechanical and chemical transitions: 

( exp {-ef) , if (i ^ j) = (1 -> 3) 
%{f) ■■= < exp ((1 - e)f) , if (i ^ j) = (3 ^ 1), 

[ l+exp [x*/] ’ 

where 0 and x] = xl additional experimental param¬ 
eters. 

For the six-state model the transitions of the forward 
cycle J- = Cl reflect experimental data. We briefly 
outline how we use the arguments of Ref. m for the 
parametrization of the four-state model shown in Fig. 
First note that transitions associated to the edges (1 


3) and (1 4) are also present in the six-state model. We 

thus use similar parameterizations. Transition (3 ^ 4) 
combines ADP release at the leading site with hydrolysis 
(and immediate P release) at the trailing one. In the six- 
state model, the same numerical value of the mechanical 
parameter x*- = 0.15 is used for both of these transitions. 
We take this as a motivation for using X 4 = xl = 0.15 
in the four-state model. Now the only undetermined 
parameters in the forward cycle are the first-order rate 
constants and The Hill-Schnakenberg condition 
a I = — / + A/i for vanishing mechanical driving / = 0 
yields one additional constraint which can be cast into 
the expression 

3 4 1 

444 A rx 

4 1 3 — -^eq • 

Finally, we take as a fit parameter that we use to ad¬ 
just our model to experimental results at the physiologi¬ 
cal parameter values: we choose it such that the six-state 
and four-state model yield the same average velocity c{d) 
for / = 0 and [ATP] = [ADP] = [P] = 1 pM. 

The parameters of the remaining transitions are ob¬ 
tained by symmetry. The exception is the first-order con¬ 
stant /^ 2 , associated with ATP release from the leading 
head. Similar to Ref. m we adjust it in order to account 
for the Hill-Schnakenberg conditions. On the dissipative 
cycle V = C 2 this constraint amounts to (J 2 = 2 A/i and 

yields 4 = (^) 4 - 

At this point, we have determined all the parameters of 
the four-state model while ensuring the physical and ther¬ 
modynamic consistency with the original six-state model. 
Fitting yields = 2.52 x 10^, which proves to be a good 
choice globally, cf. Fig.[^ All model parameters are sum¬ 
marized in Table m 
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